Impact of IDH Mutations, the 1p/19q Co-Deletion and the G-CIMP Status on Alternative Splicing in Diffuse Gliomas

By generating protein diversity, alternative splicing provides an important oncogenic pathway. Isocitrate dehydrogenase (IDH) 1 and 2 mutations and 1p/19q co-deletion have become crucial for the novel molecular classification of diffuse gliomas, which also incorporates DNA methylation profiling. In this study, we have carried out a bioinformatics analysis to examine the impact of the IDH mutation, as well as the 1p/19q co-deletion and the glioma CpG island methylator phenotype (G-CIMP) status on alternative splicing in a cohort of 662 diffuse gliomas from The Cancer Genome Atlas (TCGA). We identify the biological processes and molecular functions affected by alternative splicing in the various glioma subgroups and provide evidence supporting the important contribution of alternative splicing in modulating epigenetic regulation in diffuse gliomas. Targeting the genes and pathways affected by alternative splicing might provide novel therapeutic opportunities against gliomas.


Introduction
Recent advances in our understanding of the molecular alterations underlying diffuse gliomas have led to a novel classification based on molecular markers including isocitrate dehydrogenase (IDH) 1 and 2 mutations and 1p/19q co-deletion (codel) [1]. It is now recognized that IDH-wild-type (WT) gliomas essentially represent primary glioblastomas (GBM), while IDH-mutant gliomas bearing (1p/19q codel) or not bearing (1p/19q non-codel) the co-deletion of 1p/19q encompass tumors previously classified as oligodendrogliomas and astrocytomas, respectively. The subsequent incorporation of DNA methylation profiling, based on the discovery of a CpG island methylation phenotype (CIMP) in a subset of gliomas [2], led to a refinement in the subclassification of IDH-mutant 1p/19q non-codel gliomas in G-CIMP-high and G-CIMP-low subgroups (reviewed in [3,4]). Notably, all of the 1p/19q codel and GCIMP-high gliomas are mutated in IDH1 or IDH2 [2,5,6].
A recent analysis of the alternative splicing landscape of pediatric and adult highgrade glioma (HGG) has uncovered an increased splicing burden compared with that in normal brain [7]. However, how IDH mutations, the co-deletion of 1p/19q and the G-CIMP phenotype affect alternative mRNA splicing and the expression of protein isoforms in diffuse gliomas remains largely unknown.
In gliomas, IDH mutations are found primarily in IDH1 (incidence > 70%), with R132H representing > 90% of all mutations identified. Mutations in IDH2, affecting the analogous residue R172, have also been identified, but they are much rarer [8,9]. Of note, IDH mutations are also found in multiple other tumors [10]. At the mechanistic level, the IDH mutations define a neomorphic activity. Specifically, whereas WT IDH catalyzes the conversion of isocitrate into α-KG, mutant IDH converts α-KG into oncometabolite 2HG, which is an inhibitor of multiple αKG-dependent dioxygenases, including DNA/RNA demethylases, histone demethylases and proline/lysine hydroxylases [10,11]. In addition to metabolic alterations [12][13][14], the accumulation of 2HG impacts the removal of crucial epigenetic marks, leading to alterations in a variety of cellular programs affecting cellular metabolism, cancer biology and oncogenesis [10,14,15]. Notably, the TET family of DNA hydroxylases involved in DNA demethylation [16] is a major pathological target of the IDH mutations [10]. Thus, neomorphic IDH1-mutant activity leads to a glioma CpG island methylator phenotype (G-CIMP) characterized by concurrent promoter hypermethylation and the silencing of a large number of genes [2,6]. As for the IDH mutations [17], and contrasting with our knowledge of the impact of the G-CIMP phenotype on the DNA methylome and transcriptome in gliomas [2,4,18], very little is known about the impact of G-CIMP status on alternative mRNA splicing and the expression of protein isoforms.
The list of αKG-dependent dioxygenases acting on DNA/RNA and affected by the IDH mutations also contains crucial erasers of a fundamental epigenetic mark affecting RNA: the ALKBH5 and FTO dioxygenases that mediate oxidative removal of N6methyladenosine (m6A) in RNA [19,20]. This modification, which defines a dynamic RNA epigenetic code, is interpreted by specific readers that regulate a variety of biological processes, including the alternative splicing of precursor (pre)-mRNA-a process that dramatically increases the complexity of gene expression [21][22][23]. Epigenetic readers of the m6A-RNA mark involved in pre-mRNA splicing include HNRNPA2B1 [24], HNRNPC [25], YTHDC1 [26] and SRSF2 [27], whose interaction with RNA is modulated by RNA methylation. Compounding the problem, the activity of crucial splicing factors such as U2AF65 has been found to be directly regulated through protein post-translational modification (i.e., lysyl-5-hydroxylation) by the αKG-dependent lysine hydroxylase Jmjd6, with important biological consequences [28,29], suggesting an additional mechanism whereby 2HG might affect pre-mRNA splicing.
The evidence indicates that the 1p/19q co-deletion impacts alternative splicing at least in part through its effect on the splicing factor encoded by the FUBP1 (far upstream element-binding protein 1) gene, located on chromosome 1. Mutations in the remaining allele of FUBP1 are frequently encountered in oligodendrogliomas carrying the 1p/19q co-deletion [30]. A recent study of the somatic mutational landscape of splicing factor genes across 33 cancer types revealed a significant association of loss-of-function (LoF) mutations affecting the remaining allele of FUBP1 with alternative splicing and gene downregulation in oligodendrogliomas [31].
Here, we have carried out a bioinformatics analysis to examine the impact of the IDH mutation, as well as that of the 1p/19q co-deletion, the G-CIMP status and FUBP1 status, on alternative splicing and protein isoform expression in a cohort of 662 diffuse gliomas from The Cancer Genome Atlas (TCGA).

Splicing Alterations in Gliomas
To explore the impact of mutant IDH, as well as that of the 1p/19q co-deletion and the G-CIMP status, on alternative splicing in gliomas, we undertook an investigation of the transcriptome RNA sequencing (RNAseq) data of a cohort of 662 diffuse gliomas (see Materials and Methods, Section 4.1) contained in The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/, accessed on 6 September 2022)) [32]. The samples of this cohort were segregated into the following groups: IDH-wild-type (WT) gliomas, IDHmutant + 1p/19q non-codel gliomas, and IDH-mutant + 1p/19q codel gliomas. The latter group was further divided into FUBP1 LoF and FUBP1-WT subgroups [31], whereas the 1p/19q non-codel group was divided into G-CIMP-high and G-CIMP-low subgroups ( Figure 1A). Notably, we did not identify FUBP1 mutations in the IDH-WT group (see Supplementary Figure S5A). https://cancergenome.nih.gov/, accessed on 6 September 2022)) [32]. The samples of this cohort were segregated into the following groups: IDH-wild-type (WT) gliomas, IDHmutant + 1p/19q non-codel gliomas, and IDH-mutant + 1p/19q codel gliomas. The latter group was further divided into FUBP1 LoF and FUBP1-WT subgroups [31], whereas the 1p/19q non-codel group was divided into G-CIMP-high and G-CIMP-low subgroups (Figure 1A). Notably, we did not identify FUBP1 mutations in the IDH-WT group (see Supplementary Figure S5A).   We then used the percent spliced in index (PSI; see Materials and Methods, Section 4.1) to examine the efficiency of exon splicing into the transcript population of a given gene in the following comparisons: IDH-mutant vs. IDH-WT, 1p/19q codel vs. non-codel, G-CIMPhigh vs. G-CIMP-low and FUBP1 LoF vs. FUBP1-WT ( Figure 1B,C). Strikingly, more than 80% of the differential splicing events associated with the comparison of FUBP1 LoF vs. FUBP1-WT were exon skipping, which is in agreement with the documented role of FUBP1 in regulating exon splicing and modulating exon inclusion [31,[33][34][35], as well as its role in shaping the splicing landscape in low-grade gliomas [31]. However, no functional network could be identified from the list of genes affected by exon skipping for this comparison. In the other comparisons, the most frequently encountered differential splicing alterations again included exon skipping events, as well as alternate promoter (alternative first exon) and alternate terminator (alternative last exon) events, which together composed 80% of the observed events ( Figure 1C).
Next, we compared the lists of differentially spliced genes obtained above with the lists of differentially expressed genes (DEGs; see Materials and Methods, Section 4.2) which were obtained for the same comparisons ( Figure 1D). Supplementary Table S1 provides detailed lists of the alternative splicing (AS) events in all of the genes we found affected either by alternative splicing and differential gene expression (list 1), or by alternative splicing only (list 2), within the four comparisons. For the remaining part of this work, we decided to focus on those genes that were affected by alternative splicing but not by differential gene expression (Supplementary Table S1, list 2).
We first examined the chromosomal distribution of the genes with AS events for the four comparisons, using gene feature annotations from the Molecular Signatures Database (MSigDB) [36], as well as feature annotations from the NCBI gene_info file (ftp://ftp.ncbi.  Figure S1B). These observations suggest that, unlike for chromosome 19 which appears relatively unaffected by the 1p/19q co-deletion, an important portion of the genes affected by alternative splicing on chromosome 1 is present on the p arm of this chromosome. In the G-CIMP-high vs. G-CIMP comparison, where all samples have retained 1p and 19q, both chromosomes 1 and 19, as well as chromosomes 12, 16 and 17, were significantly enriched for AS, while several chromosomes displayed slightly fewer AS events than expected (Supplementary Figure S1C). Finally, in the FUBP1 LoF vs. FUBP1-WT comparison of 1p/19q codel samples, we did not observe significant differences in the genomic location of genes affected by AS compared to other genes, which was likely due to the low number of events (40 alternatively spliced genes) (Supplementary Figure S1D).
In summary, our data suggest that chromosome arm 1p, but not 19q, is significantly affected at the AS level by the 1p/19q co-deletion. Although FUBP1 is located on chromosome 1p, and while we cannot exclude a role for this factor in the splicing of a subset of the genes located on this chromosome arm, our data also suggest that the loss of AS events on this chromosome arm in the 1p/19q codel subgroup is not the result of FUBP1 loss of function. Finally, our observations of chromosome 19 suggest that the loss of one chromosome 19q arm in the 1p/19q codel subgroup triggers "compensatory" AS events on the intact arm, presumably participating in the development of these tumors.
We next examined the distribution of AS events as a function of the number of exons in the affected genes (see Materials and Methods, Section 4.3). The data reported in Supplementary Figure S2 indicate that, for all comparisons, fewer AS events than expected based on ShinyGo [37] were observed for genes composed of less than four exons, whereas they tended to be enriched for genes above four exons, with maximum enrichments observed for genes with five and six exons. Finally, we examined the distribution of AS events as a function of the number of transcripts associated with the differentially spliced genes. Overall, our analysis revealed significantly higher numbers of AS genes than predicted for genes expressing more than six transcripts (Supplementary Figure S3).
Focusing on AP and AT, we next attempted to trace the AS events identified across the four comparisons (IDH-mutant vs. IDH-WT, 1p/19q codel vs. non-codel, G-CIMP-high vs. G-CIMP-low and FUBP1 LoF vs. FUBP1-WT) to specific molecular alterations (IDH mutations, 1p/19q co-deletion and G-CIMP-high). To this end, we carried out a Venn diagram analysis of the data from Supplementary Table S1 (list 2) to identify shared as well as unique AP and AT events within each comparison (Supplementary Figure S4A). Such analysis allowed us to identify AS events associated with the IDH mutational status and which (i) remained unchanged in the subsequent subgroup ramifications, or (ii) were further deregulated in one or several subgroups, as well as (iii) AS events unaffected by the IDH mutations, but associated with one or several subgroups. Such AS events and their affected genes are represented in Supplementary Figure S4B and listed in Supplementary  Table S2. Only three genes were affected by AS events associated with the FUBP1 LoF vs. FUBP1-WT comparison, likely because of the small number of FUBP1 LoF specimens in our cohort. In total, we identified 207 and 327 genes linked to IDH-mutant-independent AS events (AP or AT) in the 1p/19q codel vs. non-codel and G-CIMP-high vs. G-CIMP-low comparisons, respectively. In contrast, 748 genes were linked to AS events associated solely with the IDH mutational status, which appear to remain fixed in the subgroup ramifications.

Biological Processes and Molecular Functions Affected by Alternative Splicing in Gliomas
To explore the effect of splicing changes on the biology of the various glioma subgroups, we next performed a gene ontology (GO) analysis using ShinyGo [37] on the genes with alternative splicing events identified in the four comparisons (see Materials and Methods, Section 4.3). When biological processes were considered, we found that terms related to intracellular transport, protein/macromolecule localization and RNA processing were enriched in the IDH-mutant vs. IDH-WT comparison ( Figure 2A). Remarkably, the 1p/19q codel vs. non-codel comparison revealed a strong enrichment of terms related to neuron morphogenesis, differentiation and development ( Figure 2B). Although we found that intracellular transport and protein/macromolecule localization were enriched in the G-CIMP-high vs. G-CIMP-low comparison, the G-CIMP status also appeared to affect unique biological processes at the splicing level. These included histone acetylation, as well as RNA splicing itself ( Figure 2C), the latter being illustrated by the presence of several splicing factors among the differentially spliced genes identified for this comparison (Table 1).   Previous studies have associated the alternative splicing of arginine-serine-rich (SR) splicing factors to nonsense-mediated mRNA decay [38]. However, this form of unproductive splicing cannot be the reason for our observations, as we analyzed solely alternatively spliced genes that are not affected at the level of gene expression. No significant enrichments were observed for the FUBP1 LoF vs. FUBP1-WT comparison.
We next considered the GO molecular functions. We identified that the differentially spliced genes were linked to significant molecular functions in each comparison. We found that mRNA-binding and transcription coregulatory activity were enriched in IDHmutant vs. IDH-WT ( Figure 3A) and in G-CIMP-high vs. G-CIMP-low group comparisons ( Figure 3C), but not in the 1p/19q codel vs. non-codel comparison ( Figure 3B). This may indicate that the presence of transcript isoforms may further impact specific signal transduction of RNAs. Another striking feature of all comparisons was the presence of terms related to the small GTPase activity and its regulation ( Figure 3A-C). Indeed, of all the molecular function terms found enriched in our comparisons, GTPase activity was the only common one ( Figure 3D). A list of the GTPase genes identified in our study is presented in Supplementary Table S3.

Modulation of Epigenetic Regulation by Alternate Splicing in Gliomas
In a previous comparison of pediatric high-grade glioma (pHGG) samples from the normal brain, Siddaway et al. reported that chromatin modifiers were more affected by splicing than expression changes in pHGG, implying that AS is a crucial modulator of epigenetic regulation in these gliomas [7]. To examine this issue in the diffuse glioma cohort, we compared the lists of differentially spliced and expressed genes identified in our various comparisons with that of the Database of Epigenetic Modifiers (dbEM) [39]. When the IDH-mutant vs. -WT comparison was considered, we found roughly equal numbers of epigenetic modifiers in the genes affected by splicing only (N = 38) or by expression change only (N = 35), with another 11 epigenetic modifiers affected at both levels ( Figure 4A). These observations underline the importance of AS in the epigenetic changes brought about by the IDH mutations. Notably, differential splicing and gene expression also contributed equally to the deregulation of epigenetic modifiers in the 1p/19q codel vs. non-codel and G-CIMP-high vs. G-CIMP-low comparisons, whereas only one epigenetic modifier was affected in the FUBP1 LoF vs. FUBP1-WT comparison ( Figure 4A).
We next examined which epigenetic modifiers were affected by differential splicing in our comparisons (Table 2).
We found that DNMT3A, INO80E, MTA1 and PRMT2 were common to the IDH-mutant vs. -WT, 1p/19q codel vs. non-codel and G-CIMP-high vs. G-CIMP-low comparisons (Table 2). Interestingly, DNA methyltransferase DNMT3A and two members of the methyl-CpG-binding domain (MBD) protein family, MBD1 and MBD2, were among the epigenetic modifiers affected by differential splicing in several comparisons, including the G-CIMPhigh vs. G-CIMP-low comparison. Specifically, in the case of MBD1, we found that the splicing index of exon 10 was increased in the IDH-mutant glioma (mean PSI: 0.31) compared to IDH-WT (mean PSI: 0.13), indicating the higher expression of an MBD1 isoform containing this exon in the IDH-mutant glioma. Conversely, a small decrease in the expression of an MBD1 isoform containing exon 12 was predicted in the IDH-mutant glioma ( Figure 4B). MBD1 splicing variants were also predicted for the G-CIMP-high vs. G-CIMP-low comparison where the percentage of the MBD1 isoforms containing the exons flanked by exon 17 and 18.5, i.e., 18.1-18.4, was decreased from 93% to 87% ( Figure 4B). For MBD2, which encodes two isoforms using exon 3.2 and exon 8 as terminators, respectively, we found that usage of exon 8 as the terminator was higher in the 1p/19q codel subgroup (mean PSI: 0.96) than in the non-codel subgroup (mean PSI: 0.92), while the opposite was seen for exon 3.2 (mean PSI: 0.04 vs. 0.08). A similar pattern was observed in the G-CIMPhigh vs. G-CIMP comparison, where usage of exon 8 as the terminator was higher in the G-CIMP-high subgroup than in the G-CIMP-low subgroup.  We next examined which epigenetic modifiers were affected by differential splicing in our comparisons ( Table 2).

Glioma Driver Alterations and Alternate Splicing
Having considered the impact of IDH mutations on splicing, we next examined if other mutations that have been implicated as glioma drivers or as clinically relevant markers for subgroup classification were associated with specific alternatively spliced genes. For this study, we examined EGFR, PTEN, CIC, TP53 and ATRX, as well as mutations affecting the TERT promoter [30,[41][42][43], using the cBioPortal for Cancer Genomics (see Materials and Methods, Section 4.4). The molecular alterations affecting these genes in the 662 specimens of our cohort are illustrated in Supplementary Figure S5A. Notably, EGFR alterations (mutations and amplifications, considered together in this study) and PTEN alterations (mutations and deep deletions) affected almost exclusively IDH-WT gliomas. On the other hand, mutations affecting CIC and FUBP1 were restricted to the IDH-mutant 1p/19q codel subgroup. ATRX mutations were almost exclusively present in the IDH-mutant 1p/19q non-codel subgroup, and this subgroup also contained the majority of the TP53 mutations. TERT promoter mutations were distributed across all subgroups. We next identified the AS genes associated with EGFR, PTEN, CIC, TP53, ATRX and TERT. Because alterations affecting these genes were strongly confined to specific subgroups (Supplementary Figure S5A), and to prevent interference from other subgroups, we restricted each analysis to the subgroup(s) hosting the gene alteration under consideration. The identified AS genes are represented in Supplementary Figure S5B and listed in Supplementary Table S4. Very few AS genes were associated with alterations in EGFR, PTEN, CIC and ATRX. With the exception of TP53, the GO analysis using DAVID did not reveal any significant enrichment in the biological processes and molecular functions for the identified AS genes associated with the driver genes considered in this work. However, when the mutational status of TP53 was considered within the IDH-mutant 1p/19q non-codel subgroup, the GO analysis of the AS genes revealed the enrichment of molecular function terms related to small GTPase-binding (GO:0031267; p-value 1.4 × 10 −8 , FDR 2.7 × 10 −6 ) as well as GTPase-activator activity (GO:0005096; p-value 7.28 × 10 −4 , FDR 0.04).
Finally, the most frequently encountered differential splicing alterations associated with TP53 mutations in the IDH-mutant 1p/19q non-codel subgroup included exon skipping events (53.1%), whereas alternate promoter (37.1%) and exon skipping (35.3%) predominated in the TERT-mutant vs. -WT comparison. Notably, alternate promoter events represented almost two thirds of the AS events in the EGFR-mutant vs. -WT comparison, where no alternate terminator event was detected (Supplementary Figure S5C).

Discussion
The contribution of alternative splicing to cancer development and biology has received increasing attention in recent years [44,45]. Previous studies have identified alternative splicing as an important oncogenic pathway in high-grade gliomas [7]. In addition, mRNA splicing profiling has been considered as an approach to identify prognostic predictors and potential therapeutic targets for GBM [46]. In this work, we have investigated the impact of IDH mutations, the co-deletion of 1p/19q and the G-CIMP phenotype on alternative mRNA splicing and the expression of protein isoforms in diffuse gliomas.
Our analysis reveals long lists of genes affected at the splicing level by the IDH mutations, the co-deletion of 1p/19q and the G-CIMP phenotype. They also suggest that some of the AS events associated with the IDH mutations remain fixed in the subsequent ramifications afforded by the 1p/19q codel/non-codel and G-CIMP-high/low phenotypes, while others were specifically associated with these phenotypes. Such analyses should prove a useful resource for readers who want to examine splicing alterations in specific genes and/or clinical subgroups. Notably, several genes and gene families in this list are involved in normal neural/glial function and/or glioma proliferation, adding weight to the notion that pre-mRNA splicing alterations have a significant impact on the biology and behavior of gliomas.
Strikingly, GTPase activity and its regulation are crucial molecular functions affected by alternative splicing in diffuse gliomas. GTPases are small GTP-binding proteins that operate as molecular switches in cytoskeletal dynamics, cell communication, intracellular trafficking and cell migration [47]. GTPases exist in either an active GTP-bound state, during which the GTPase can interact with downstream effectors, or an inactive GDPbound state. Cycling between these two states is regulated by guanine nucleotide exchange factors, GTPase-activating proteins and guanine nucleotide dissociation inhibitors [48,49]. Small GTPase genes play a crucial role in brain development [50,51]. Small GTPases have been found to be involved in GBM motility, invasion and progression [52,53]. In addition, the evidence indicates that Rho family GTPases modulate DNA repair and resistance to ionizing radiation in WT p53 GBM [54]. GTPases are considered as therapeutic targets in cancer, including gliomas [55][56][57]. The small GTPase family is encoded by dozens of genes displaying complex splicing dynamics [58]. Our data suggest that the repertoire of GTPases and their regulators is subjected to intensive modulation by alternative splicing in key subgroups of diffuse gliomas. How these alternative splicing events contribute to the development, biology, treatment response and prognosis of the various glioma subgroups examined in this study remains to be investigated. In this regard, our data suggest that TP53 mutations are associated with AS events affecting GTPase-binding and regulation factors in the IDH-mutant 1p/19q non-codel subgroup. Interestingly, altered RNA splicing by the mutant p53, through the modulation of the RNA-binding protein hnRNPK, has been shown to impact GTPase-activating proteins, leading to the activation of oncogenic RAS signaling in pancreatic cancer [59]. Whether the biological processes identified in our analysis reflect a similar role for p53 in splicing regulation and oncogenesis in specific glioma subgroups remains unknown.
Alternative splicing acts as a crucial modulator of epigenetic regulation in pHGG [7]. We found evidence suggesting that this role can be extended to all diffuse gliomas. Remarkably, among the epigenetic modifiers affected by differential splicing in our comparisons were DNA methyltransferases and epigenetic readers of DNA methylation. Although the functional impact of the alterations conferred by the differential splicing of these factors remains to be determined, these observations suggest that the IDH mutations elicit an integrated reprogramming that affects not only DNA methylation patterns, but also the expression of variant epigenetic readers of DNA methylation.
Several strategies targeting splicing mechanisms, splicing factors or alternative splicing isoforms are being considered in relation to cancer therapy [44,45]. Our analyses reveal that a great number of genes involved in glioma-relevant biological pathways are affected by alternative splicing in diffuse gliomas. Targeting these genes might provide novel therapeutic opportunities against gliomas.

Differential Splicing Analysis
We downloaded the PSI data pertaining to the diffuse glioma cohort (670 patients with diffuse glioma, including 515 low-grade glioma (LGG) and 155 glioblastoma (GBM), from the TCGA SpliceSeq database (https://bioinformatics.mdanderson.org/TCGASpliceSeq/ PSIdownload.jsp, accessed on 6 September 2022) [60] as well as that of molecular subtypes (clinical subgroup, i.e., IDH status, 1p/19q co-deletion and Glioma CpG island methylator phenotype) of the TCGA samples using the TCGAquery_subtype function of the TCGAbiolinks R package (10.1371/journal.pcbi.1006701) and the FUBP1 mutation data from the TCGA MC3 working group [61] using the getMC3MAF function of the TCGAbiolinks R package. Then, 8 samples with IDH information missing were removed and the data from the remaining 662 samples were used for further analysis. We performed differential expression analysis according to the methods described by Seiler et al. [31]. Briefly, we filtered out the splicing sites which had a 0 PSI in more than half of the samples from both groups in the comparison. We converted each PSI to log odds (log(PSI/(1-PSI))), and then fit a linear model for each splicing site and computed the moderate t-statistic to define differential expression using the limma R package [62]. We considered the splicing sites with adjusted p-values (FDR: false discovery rate) < 0.05 and log2 fold changes (FC) > 0.5 or <−0.5 as significantly regulated for all types of splicing events.

Differential Gene Expression Analysis
We downloaded the RNA-seq raw counts of the GBM (N = 173) and LGG (N = 529) samples from the TCGA database (https://portal.gdc.cancer.gov/, accessed on 4 February 2018). We kept the data from 656 primary tumor samples, which also had PSI data, for further analysis. We filtered out the genes detected with less than 10 reads in more than half of the samples. We then used the functions DESeq and results from the DESeq2 R package [63] to perform size factor estimation, dispersion estimation and negative binomial GLM fitting sequentially for differentially expressed gene identification. We considered the genes with an FDR < 0.05 and a log2 FC change > 0.5 or <−0.5 to be differentially expressed.

Over-and Under-Representation Analyses and Pathway Enrichment
In general, the over-representation analysis was conducted using ShinyGO 0.77 (http://bioinformatics.sdstate.edu/go, accessed on 30 March 2023/) [37], which allowed for direct visual representation, or DAVID (Database for Annotation, Visualization and Integrated Discovery) [64,65]. For diagrams we used Pathview [66] with KEGG database pathway annotation [67]. A cut-off of an FDR < 0.01 was applied for significance. For Supplementary Figure S1, the over-or under-representation of individual chromosomes was analyzed using the hypergeometric test. Gene locations were obtained with feature annotations from the NCBI database (ftp://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/ Mammalia/Homo_sapiens.gene_info.gz, accessed on 23 May 2023). Analysis was performed in R software (R.4.2.3) using the phyper function, where the lower.tail parameter was set to TRUE for under-representation and FALSE for lower over-representation. p-values were adjusted for multiple testing errors using the Benjamini-Hochberg procedure [68]. All genes present in the database were considered as the background data.